MODULE datamm
  USE GLOBVAR
  USE STATISTICS
  !USE INTEGRATION
  IMPLICIT NONE
  
CONTAINS
  
  SUBROUTINE datam(momin)
    IMPLICIT NONE
    REAL(8), INTENT(INOUT):: momin(:) !initial guess
    REAL(8) :: a(N_ind,n_testl)
    REAL(8) :: inter(n_ind,n_levell,20,3),interl(n_ind,n_levell,3)
    REAL(8),ALLOCATABLE:: AVSE(:),cor(:,:)
    INTEGER :: du(n_levell),ind(N_ind,n_levell),mmm,nm,mn 
    INTEGER :: i,j,k,t,m,n,l,indexl(n_testl,2), artl(n_ind,11,20),con(n_mom)
    
    !!!!! start to calculate data moments!!!!!
    
    !!! language tasks!!!!
    indexl=0
    DO i=1,N_ind
      DO j=1,N_testl
         IF (l_test(i,j)>-90) THEN
             indexl(j,1)=indexl(j,1)+1
             IF (l_test(i,j)>0.5) indexl(j,2)=indexl(j,2)+1
         END IF
         
      END DO
    END DO
    
    DO i=1,N_testl
        IF (indexl(i,1)>0.5) THEN
            momin(i)=DBLE(indexl(i,2))/DBLE(indexl(i,1)) ! PASSING RATE FOR EACH LANGUAGE TASK
            con(i)=indexl(i,1)
        ELSE 
            momin(i)=-99
            con(i)=-99
        END IF
    END DO
    
    !!!! generate artl variables !!!!
    artl=-1
    DO i=1,N_ind
        DO j=1,N_testl
            IF (l_test(i,j)>-90) THEN
                k=levell(j)
                m=repl(i,j)
                artl(i,k,m)=l_test(i,j)
                inter(i,k,m,1:3)=x_testl(i,j,1:3)
            END IF
        END DO
    END DO
    
    !!!! passing rate at each level!!!!

    DO j=N_testl+1,N_testl+N_levell
       mn=0
       nm=0

       DO i=1,N_ind
           m=0
           n=0 
           DO k=1,20
               IF (artl(i,j-N_testl,k)>-0.5) THEN
                   m=m+1
                   mn=mn+1
                   IF (artl(i,j-N_testl,k)==1) THEN
                       n=n+1
                       nm=nm+1
                   END IF
                   
               END IF
           END DO 
           IF (m>0) THEN
               ave(i,j-N_testl)=DBLE(n)/DBLE(m)
           ELSE
               ave(i,j-N_testl)=-99
           END IF
       END DO
       IF (mn>0) THEN
           momin(j)=DBLE(nm)/DBLE(mn)
           con(j)=m
       ELSE
           momin(j)=-99
           con(j)=-99
       END IF
       
    END DO
    
    


    
    !!!! correlation across levels !!!!
    
    DO j=2,N_levell-1
      DO t=1,5
        DO k=1,5
          m=0
          n=0 
          DO i=1,N_ind 
             IF (artl(i,j,t)==1 .and. artl(i,j+1,k)>-0.5) THEN    ! the probability of passing the kth task at level j+1 conditional on passing the tth task at level j
                 m=m+1
                 IF (artl(i,j,t)==1 .and. artl(i,j+1,k)==1) n=n+1
             END IF
          END DO
          
          IF (m>0) THEN
             momin(N_testl+N_levell+(j-2)*25+(t-1)*5+k)=DBLE(n)/DBLE(m)
             con(N_testl+N_levell+(j-2)*25+(t-1)*5+k)=m
          ELSE
             momin(N_testl+N_levell+(j-2)*25+(t-1)*5+k)=-99
             con(N_testl+N_levell+(j-2)*25+(t-1)*5+k)=-99
          END IF
        END DO
      END DO
    END DO
    
    DO j=2,N_levell
        l=0
      DO t=1,4
        DO k=t+1,5
          m=0
          n=0
          l=l+1
          DO i=1,N_ind 
             IF (artl(i,j,t)==1 .and. artl(i,j,k)>-0.5) THEN    ! the probability of passing the kth task at level j+1 conditional on passing the tth task at level j
                 m=m+1
                 IF (artl(i,j,t)==1 .and. artl(i,j,k)==1) n=n+1
             END IF
          END DO
          
          IF (m>0) THEN
             momin(25*(N_levell-2)+N_testl+N_levell+(j-2)*10+l)=DBLE(n)/DBLE(m)
             con(25*(N_levell-2)+N_testl+N_levell+(j-2)*10+l)=m
          ELSE
             momin(25*(N_levell-2)+N_testl+N_levell+(j-2)*10+l)=-99
             con(25*(N_levell-2)+N_testl+N_levell+(j-2)*10+l)=-99
          END IF
        END DO
      END DO
    END DO
    
    DO j=2,N_levell-2
      DO t=1,5
        DO k=1,5
          m=0
          n=0 
          DO i=1,N_ind 
             IF (artl(i,j,t)==1 .and. artl(i,j+2,k)>-0.5) THEN    ! the probability of passing the kth task at level j+1 conditional on passing the tth task at level j
                 m=m+1
                 IF (artl(i,j,t)==1 .and. artl(i,j+2,k)==1) n=n+1
             END IF
          END DO
          
          IF (m>0) THEN
             momin(25*(N_levell-2)+N_testl+N_levell+(N_levell-1)*10+(j-2)*25+(t-1)*5+k)=DBLE(n)/DBLE(m)
             con(25*(N_levell-2)+N_testl+N_levell+(N_levell-1)*10+(j-2)*25+(t-1)*5+k)=m
          ELSE
             momin(25*(N_levell-2)+N_testl+N_levell+(N_levell-1)*10+(j-2)*25+(t-1)*5+k)=-99
             con(25*(N_levell-2)+N_testl+N_levell+(N_levell-1)*10+(j-2)*25+(t-1)*5+k)=-99
          END IF
        END DO
      END DO
    END DO
    
    indexl=0
    DO i=1,N_ind
      DO j=1,N_testl
          a(i,j)=firstl(i,2)-monthl(j)-DBLE(timel(j))/4.0d0
         IF (l_test(i,j)>-90 .and. a(i,j)>-1.25d0) THEN
             indexl(j,1)=indexl(j,1)+1
             IF (l_test(i,j)>0.5) indexl(j,2)=indexl(j,2)+1
         END IF
         
      END DO
    END DO
    
    DO i=1,N_testl
        IF (indexl(i,1)>0.5) THEN
            momin(25*(2*N_levell-5)+N_testl+N_levell+(N_levell-1)*10+i)=DBLE(indexl(i,2))/DBLE(indexl(i,1)) ! PASSING RATE FOR EACH LANGUAGE TASK
            con(25*(2*N_levell-5)+N_testl+N_levell+(N_levell-1)*10+i)=indexl(i,1)
        ELSE 
            momin(25*(2*N_levell-5)+N_testl+N_levell+(N_levell-1)*10+i)=-99
            con(25*(2*N_levell-5)+N_testl+N_levell+(N_levell-1)*10+i)=-99
        END IF
    END DO
    
    indexl=0
    DO i=1,N_ind
      DO j=1,N_testl
          a(i,j)=DBLE(firstl(i,2)-monthl(j)-DBLE(timel(j))/4.0d0)

          
         IF (l_test(i,j)>-90 .and. a(i,j)<-1.25d0) THEN
             indexl(j,1)=indexl(j,1)+1
             IF (l_test(i,j)>0.5) indexl(j,2)=indexl(j,2)+1
         END IF
         
      END DO
    END DO


    
    DO i=1,N_testl
        IF (indexl(i,1)>0.5) THEN
            momin(25*(2*N_levell-5)+2*N_testl+N_levell+(N_levell-1)*10+i)=DBLE(indexl(i,2))/DBLE(indexl(i,1)) ! PASSING RATE FOR EACH LANGUAGE TASK
            con(25*(2*N_levell-5)+2*N_testl+N_levell+(N_levell-1)*10+i)=indexl(i,1)
        ELSE 
            momin(25*(2*N_levell-5)+2*N_testl+N_levell+(N_levell-1)*10+i)=-99
            con(25*(2*N_levell-5)+2*N_testl+N_levell+(N_levell-1)*10+i)=-99
        END IF
    END DO
    
    l=0
    DO j=2,N_levell
      DO t=1,5
          m=0
          n=0
          l=l+1
          DO i=1,N_ind 
             IF (artl(i,j,t)>-0.5) THEN   
                 m=m+1
                 IF (artl(i,j,t)==1) n=n+1
             END IF
          END DO
          
          IF (m>0) THEN
             momin(25*(N_levell-2)+3*N_testl+N_levell+(N_levell-1)*10+25*(N_levell-3)+l)=DBLE(n)/DBLE(m)
             con(25*(N_levell-2)+3*N_testl+N_levell+(N_levell-1)*10+25*(N_levell-3)+l)=m
          ELSE
             momin(25*(N_levell-2)+3*N_testl+N_levell+(N_levell-1)*10+25*(N_levell-3)+l)=-99
             con(25*(N_levell-2)+3*N_testl+N_levell+(N_levell-1)*10+25*(N_levell-3)+l)=-99
          END IF
      END DO
    END DO
    
    mmm=25*(N_levell-2)+3*N_testl+N_levell+(N_levell-1)*10+25*(N_levell-3)+5*(N_levell-1)
    
    
    DUL=-99
    ind=-99
    DO i=1,N_ind
        DO k=2,N_levell
            DO j=1,20
                IF (artl(i,k,j)==1 .and. j==1) THEN
                    dul(i,k)=1
                    ind(i,k)=1
                ELSE IF (artl(i,k,j)==0 .and. j==1) THEN
                    dul(i,k)=j
                    ind(i,k)=0
                ELSE IF (j>1 ) THEN
                    IF (artl(i,k,j-1)==0 .and. artl(i,k,j)==0 .and. ind(i,k)==0) THEN
                      dul(i,k)=j
                       ind(i,k)=0
                       IF (j<20 .and. artl(i,k,j+1)==-1) THEN
                           dul(i,k)=j+1
                           ind(i,k)=1
                       END IF                      
                           
                    ELSE IF (artl(i,k,j-1)==0 .and. artl(i,k,j)==0 .and. ind(i,k)==1) THEN
                        dul(i,k)=dul(i,k)
                        ind(i,k)=1
                    ELSE IF (artl(i,k,j-1)==0 .and. artl(i,k,j)==1 .and. ind(i,k)==0) THEN
                        dul(i,k)=j
                        ind(i,k)=1
                    ELSE IF (artl(i,k,j-1)==0 .and. artl(i,k,j)==1 .and. ind(i,k)==1) THEN
                        dul(i,k)=dul(i,k)
                        ind(i,k)=1
                    END IF
                END IF
            END DO           
                
        END DO
    END DO
    
    DO k=1,n_levell
        m=0
        n=0
        DO i=1,n_ind
            IF (dul(i,k)>-90) THEN
               m=m+1
               n=n+dul(i,k)
            END IF
        END DO
        IF (m>0) THEN
          momin(mmm+k)=DBLE(n)/DBLE(m)
          con(mmm+k)=1 !m
        ELSE 
           momin(mmm+k)=-99
           con(mmm+k)=-99
        END IF
    END DO
    
    mmm=25*(N_levell-2)+3*N_testl+N_levell+(N_levell-1)*10+25*(N_levell-3)+5*(N_levell-1)+n_levell
    
    du=0
    DO k=1,n_levell
        DO i=1,N_ind
            IF (dul(i,k)>-90) THEN
                du(k)=du(k)+1
            END IF
        END DO
    END DO
    DO k=1,n_levell
        allocate(avse(du(k)))
        m=0
        DO i=1,n_ind
            IF (dul(i,k)>-90) THEN
                m=m+1
                avse(m)=dul(i,k)
            END IF
        END DO
        momin(mmm+k)=STDEV_V(AVSE)
        CON(MMM+K)=1
        DEALLOCATE(avse)
    END DO
    mmm=25*(N_levell-2)+3*N_testl+N_levell+(N_levell-1)*10+25*(N_levell-3)+5*(N_levell-1)+2*n_levell
    

    IND=0
    interl=0.0d0
    DO k=2,N_levell
        DO i=1,N_ind
          DO j=1,20
              IF (inter(i,k,j,1)>-90) THEN
                  ind(i,k)=ind(i,k)+1
                  interl(i,k,1)=interl(i,k,1)+inter(i,k,j,1)
                  interl(i,k,2)=interl(i,k,2)+inter(i,k,j,2)
                  interl(i,k,3)=interl(i,k,3)+inter(i,k,j,3)
              END IF
          END DO
          IF (ind(i,k)>0) THEN
               interl(i,k,1)=DBLE(interl(i,k,1))/DBLE(ind(i,k))
               interl(i,k,2)=DBLE(interl(i,k,2))/DBLE(ind(i,k))
               interl(i,k,3)=DBLE(interl(i,k,3))/DBLE(ind(i,k))
          END IF
        END DO
    END DO
    
    
    du=0
    DO k=1,n_levell
        DO i=1,N_ind
            IF (ind(i,k)>0 .AND. dul(i,k)>-90) THEN
            du(k)=du(k)+1          
            END IF
        END DO
    END DO
    
    DO k=1,n_levell
        allocate(cor(du(k),4))
        m=0
        DO i=1,n_ind
             IF (ind(i,k)>0 .AND. dul(i,k)>-90) THEN
                m=m+1
                cor(m,1)=interl(i,k,1)
                cor(m,2)=interl(i,k,2)
                cor(m,3)=interl(i,k,3)
                cor(m,4)=dul(i,k)
             END IF
        END DO
        IF (k>1) THEN
        momin(mmm+3*(k-1)+1)=Correlation(cor(:,1),cor(:,4))
        momin(mmm+3*(k-1)+2)=Correlation(cor(:,2),cor(:,4))
        momin(mmm+3*(k-1)+3)=Correlation(cor(:,3),cor(:,4))
        CON(MMM+3*(k-1)+1)=1 !du(k)
        CON(MMM+3*(k-1)+2)=1 !du(k)
        CON(MMM+3*(k-1)+3)=1 !du(k)
        ELSE 
            momin(mmm+3*(k-1)+1)=-99
            momin(mmm+3*(k-1)+2)=-99
            momin(mmm+3*(k-1)+3)=-99
            CON(MMM+3*(k-1)+1)=1
            CON(MMM+3*(k-1)+2)=1
            CON(MMM+3*(k-1)+3)=1
        END IF
        DEALLOCATE(cor)
    END DO
    mmm=25*(N_levell-2)+3*N_testl+N_levell+(N_levell-1)*10+25*(N_levell-3)+5*(N_levell-1)+5*n_levell
            
    
                
    
        
      
        
    
    countnum=con
      
          OPEN(12,file="data_initial.out")
             DO i=1,N_mom
              WRITE(12,'(2000F32.12)') momin(i)
             END DO
           CLOSE(12)    
    
    
    
    
  END SUBROUTINE datam
  
END MODULE datamm